Project

rm(list=ls())
suppressPackageStartupMessages({
library(rjags)
library(bayesplot)
library(jagsUI)
library(glue)
library(loo)
library(knitr)
library(kableExtra)
  }
)

Bayesian Learning & Monte-Carlo Simulation final project

Introduction

Setting up environment, importing the dataset and visualising a graph of the two time series we will have to fit.

rm(list=ls())

# If set to TRUE higher order models will also be trained
intense = TRUE

# New Family Houses Sold: United States
# Source: https://fred.stlouisfed.org/series/HSN1F
data      <- read.csv("gdp_inflation.csv",header=T)
data$DATE <- as.Date(data$DATE)

# Plotting time series for visualisation purposes
plot(
  data$DATE[1:305],
  data$GDP_PC1[1:305],
  type="l",
  xlab="",
  ylab="",
  main="Percent change YoY (quarterly)"
)
lines(
  data$DATE[1:305],
  data$CPIAUCSL_PC1[1:305],
  type="l",
  col="red"
)

legend(x="bottom",legend=c("GDP","CPI"),col=c("black","red"),lty=1)

At this point, we setup our environment by converting the values in both time series into numericals, and then splitting them into training and testing. Finally, we define some useful variables which we will be using through out our script.

# Dataset extraction
x_axis      <- data$DATE[1:305]
series_GDP  <- as.numeric(data$GDP_PC1[1:305])
series_CPI  <- as.numeric(data$CPIAUCSL_PC1[1:305])

# Parameter setup (RMK: we assume both time series of same length)
TRAIN_PERC  <- 0.15
N_tot       <- length(series_GDP)
N_test      <- floor(TRAIN_PERC * N_tot)
N_train     <- N_tot - N_test

To further inspect our data we can also plot the respective auto-correlation functions:

acf(series_GDP, lag.max=100)

acf(series_CPI, lag.max=100)

The auto-correlation functions of the two proceses do not clearly match any given noticeable pattern, hence a more sophisticated model selection will be needed, where multiple models will be tested and evaluated one against the other.

Utility functions

For our convenience, we define a few utility functions which will be useful later. The first one automatises the interactions with JAGS, the second one will draw plots of both the whole time-series and its out-of-sample portion only, the third one will specifically plot a time-series along with the anomalies which have been detected through out its history, and so on.

fit_JAGS <- function(
  model_string,   # JAGS model string
  data_list,      # List of input values for JAGS
  to_save,        # Vector with parameter names which should be saved
  n_iter=5000,    # Number of iterations. Should be greater than n_burnin
  n_adapt=1000,   # JAGS n.adapt internal parameter
  n_chain=1,      # Number of chains to run
  n_burnin=1000,  # Burn-in value
  n_thin=1,       # Thinning value
  verbose=FALSE   # Displays JAGS logging / output (turn on for debugging)
) {
  
  output <- jags(
    model.file=textConnection(model_string),
    data=data_list,
    parameters.to.save=to_save,
    n.adapt=n_adapt,
    n.iter=n_iter,
    n.chains=n_chain,
    n.burnin=n_burnin,
    n.thin=n_thin,
    verbose=verbose
  )
  
  return(output)
}

visualise <- function(
  x_axis,           # Values to be displayed on the x-axis of the plots
  original_series,  # Original time series with all samples
  output,           # Output from JAGS, with in-sample and out-of-sample predictions
  title="",         # Title to display on the graphs
  keep=5            # Training sampled to keep in the test graph
) {
  # Length of total series
  N <- length(x_axis)
  
  
  # Plotting entire graph (training-set + predictions)
  # In-sample predicted values with their 2.5%-97.5% quantiles
  yp_in_mean  <- output$mean$yp_in
  yp_in_q1    <- output$q2.5$yp_in
  yp_in_q2    <- output$q97.5$yp_in
  in_length   <- length(yp_in_mean)
  
  # Out-of-sample predicted values with their 2.5%-97.5% quantiles and 25%-75% quantiles
  yp_out_mean <- output$mean$yp_out
  yp_out_q1   <- output$q2.5$yp_out
  yp_out_q2   <- output$q97.5$yp_out
  yp_out_q3   <- output$q25$yp_out
  yp_out_q4   <- output$q75$yp_out
  
  # Removing potentially non-finite values from quantiles
  yp_in_q1[is.infinite(yp_in_q1)]   <- NA
  yp_in_q2[is.infinite(yp_in_q2)]   <- NA
  yp_out_q1[is.infinite(yp_out_q1)] <- NA
  yp_out_q2[is.infinite(yp_out_q2)] <- NA
  yp_out_q3[is.infinite(yp_out_q3)] <- NA
  yp_out_q4[is.infinite(yp_out_q4)] <- NA
  
  # Computing extreme values for first graph
  min_y <- min(original_series, yp_in_q1, yp_out_q1)
  max_y <- max(original_series, yp_in_q2, yp_out_q2)
  
  # Plotting REAL time series (as points only), both training set and test-set
  plot(
    x_axis,
    original_series,
    main=title, col="red",
    xlab="time", ylab="price",
    ylim=c(min_y, max_y)
  )
  lines(
    x_axis,
    original_series,
    col="red"
  )
  
  # Plotting a vertical line signalling end of training set and horizontal
  # Showing the computer mean
  abline(v=x_axis[in_length], col="magenta", lwd=1, lty=2)
  abline(h=output$mean$m0,    col="green",   lwd=1, lty=2)
  
  # Plotting in-sample predictions as well as 2.5%-97.5% quantiles
  points(x_axis[1:in_length], yp_in_mean, col="blue", pch="*")
  lines(x_axis[1:in_length],  yp_in_q1,   col="blue", lwd=1.5)
  lines(x_axis[1:in_length],  yp_in_q2,   col="blue", lwd=1.5)
  
  # Plotting out-of-sample predicted values as well as 2.5%-97.5%
  # quantiles and 25%-75% quantiles
  points(x_axis[(in_length+1):N], yp_out_mean, col="orange", pch="*")
  lines(x_axis[(in_length+1):N],  yp_out_q1,   col="orange", lwd=1.5)
  lines(x_axis[(in_length+1):N],  yp_out_q2,   col="orange", lwd=1.5)
  lines(x_axis[(in_length+1):N],  yp_out_q3,   col="orange", lwd=1.5, lty=2)
  lines(x_axis[(in_length+1):N],  yp_out_q4,   col="orange", lwd=1.5, lty=2)
  

  # Plotting onyl out-of-sample portion of graph (with last 5 training samples)
  # Re-computing extreme values for second graph
  min_y <- min(original_series[(in_length+1 - keep):N], yp_out_q1)
  max_y <- max(original_series[(in_length+1 - keep):N], yp_out_q2)
  
  # Plotting REAL time series (as points only), both training set and test-set
  plot(
    x_axis[(in_length+1 - keep):N],
    original_series[(in_length+1 - keep):N],
    main=title, col="red",
    xlab="time", ylab="price",
    ylim=c(min_y, max_y)
  )
  lines(
    x_axis[(in_length+1-keep):N],
    original_series[(in_length+1-keep):N],
    col="red"
  )
  
  # Plotting a vertical line signalling end of training set
  abline(v=x_axis[in_length+1], col="magenta", lwd=1, lty=2)
  abline(h=output$mean$m0,      col="green",   lwd=1, lty=2)
  
  # Plotting out-of-sample predicted values as well as 2.5%-97.5%
  # quantiles and 25%-75% quantiles
  points(x_axis[(in_length+1):N], yp_out_mean, col="orange", pch="*")
  lines(x_axis[(in_length+1):N],  yp_out_q1,   col="orange", lwd=1.5)
  lines(x_axis[(in_length+1):N],  yp_out_q2,   col="orange", lwd=1.5)
  lines(x_axis[(in_length+1):N],  yp_out_q3,   col="orange", lwd=1.5, lty=2)
  lines(x_axis[(in_length+1):N],  yp_out_q4,   col="orange", lwd=1.5, lty=2)
}

# Function in charge of visualising VAR models (requires custom rendering).
# This acts as a wrapper function to the standard visualise(...) function, defined above.
visualise_VAR = function(
  x_axis,               # Values to be displayed on the x-axis of the plots
  original_series_gdp,  # Original GDP time series with all samples
  original_series_infl, # Original INFLATION time series with all samples
  output,               # Output from VAR JAGS, with in-sample and out-of-sample predictions
  title="VAR",          # Title to display on the graphs
  keep=5                # Training sampled to keep in the test graph
) {
  
  # Marginalise GDP
  output_gdp = list()
  output_gdp$mean$yp_in   = output$mean$yp_in_gdp
  output_gdp$q2.5$yp_in   = output$q2.5$yp_in_gdp
  output_gdp$q97.5$yp_in  = output$q97.5$yp_in_gdp
  output_gdp$mean$yp_out  = output$mean$yp_out_gdp
  output_gdp$q2.5$yp_out  = output$q2.5$yp_out_gdp
  output_gdp$q97.5$yp_out = output$q97.5$yp_out_gdp
  output_gdp$q25$yp_out   = output$q25$yp_out_gdp
  output_gdp$q75$yp_out   = output$q75$yp_out_gdp
  output_gdp$mean$m0      = output$mean$m_gdp
  
  # Marginalise CPI
  output_infl = list()
  output_infl$mean$yp_in   = output$mean$yp_in_infl
  output_infl$q2.5$yp_in   = output$q2.5$yp_in_infl
  output_infl$q97.5$yp_in  = output$q97.5$yp_in_infl
  output_infl$mean$yp_out  = output$mean$yp_out_infl
  output_infl$q2.5$yp_out  = output$q2.5$yp_out_infl
  output_infl$q97.5$yp_out = output$q97.5$yp_out_infl
  output_infl$q25$yp_out   = output$q25$yp_out_infl
  output_infl$q75$yp_out   = output$q75$yp_out_infl
  output_infl$mean$m0      = output$mean$m_infl
  
  # Call visualisation on marginalisations
  visualise(
    x_axis,
    original_series_gdp,
    output_gdp,
    title=glue("GDP series, {title}"),
    keep
  )
  visualise(
    x_axis,
    original_series_infl,
    output_infl,
    title=glue("CPI series, {title}"),
    keep
  )
}

# Function in charge of visualising the anomaly points on top of its corresponding time series
visualise_anomalies <- function(
  x_axis,             # X-axis to use for the plot
  original_series,    # Original time series to display
  high_var_indices,   # List of indices corresponding to high variance points
  med_var_indices,    # List of indices corresponding to medium variance points
  title="",           # Name of the time series used
  start_point=1       # Value after which plotting should start (in order to obtain close up graphs)
) {
  # Plot actual time series
  N <- length(x_axis)
  plot(
    x_axis[start:N],
    original_series[start:N],
    main=glue("Anomalies for {title}"),
    t="l", xlab="time", ylab="price"
  )

  # Plotting strong anomaly points
  points(x_axis[high_var_indices], original_series[high_var_indices], col="red")
  for (high_anomaly in high_var_indices) {
    abline(v=x_axis[high_anomaly], col="red", lwd=1.5)
  }
  
  # Plotting medium anomaly points
  points(x_axis[med_var_indices], original_series[med_var_indices], col="orange")
  for (medium_anomaly in med_var_indices) {
    abline(v=x_axis[medium_anomaly], col="orange", lty=3, lwd=1.5)
  }
}

# Function in charge of classifying anomaly points given their delta value after running the anomaly detection simulation
classify_anomalies <- function(
  output,             # JAGS simulation output
  N_tot,              # Total number of samples
  title="",           # Name of time series to display
  high_value_thr=0.2, # High variance threshold value used
  med_value_thr=0.60  # Medium variance threshold value used
) {
  # Obtaining precision, variance and delta
  tau     <- output$mean$tau
  sigma2  <- 1 / tau
  delta   <- output$mean$d
  
  # Classifying the values of d
  high_var_indices  <- seq(1:(N_tot-1))[delta < high_value_thr]
  med_var_indices   <- seq(1:(N_tot-1))[high_value_thr < delta & delta < med_value_thr]
  low_var_indices   <- seq(1:(N_tot-1))[delta > med_value_thr]
  
  # Logging counts
  n_high  <- length(high_var_indices)
  n_med   <- length(med_var_indices)
  print(
    glue("Found {n_high} strong and {n_med} medium outliers out {N_tot} total samples for {title}")
  )
  
  # Returning the indices
  return(list(
    high_var_indices = high_var_indices,
    med_var_indices  = med_var_indices,
    low_var_indices  = low_var_indices
  ))
}

# Function in charge of plotting the anomaly point delta values with their relative threshold lines
display_classification <- function(
  x_axis,             # X-axis to use for the plot
  delta,              # Sequence of delta values
  high_var_indices,   # List of indices corresponding to high variance points
  med_var_indices,    # List of indices corresponding to medium variance points
  low_var_indices,    # List of indices corresponding to low variance points
  title="",           # Name of the time series used
  high_value_thr=0.2, # High variance threshold value used
  med_value_thr=0.60  # Medium variance threshold value used
) {
  # Plotting the posterior values of delta and classifying them
  plot(
    x_axis[low_var_indices], delta[low_var_indices],
    xlab="time",
    ylab="Posterior delta values",
    main=glue("{title} anomaly samples"),
    xlim=c(min(x_axis), max(x_axis)),
    ylim=c(0, 1)
  )
  
  points(x_axis[high_var_indices], delta[high_var_indices], col="red")
  points(x_axis[med_var_indices], delta[med_var_indices]  , col="orange")
  abline(h=high_value_thr, col="red")
  abline(h=med_value_thr, col="orange")
}

# K is the number of effective params
compute_metrics <- function(output, series, N_train, k) {
  
  ############################################################
  ###################### IN-SAMPLE ###########################
  ############################################################
  # In-sample: DIC, WAIC, BIC
  ll_mat <- output$sims.list$loglik
  DIC_val <- compute_dic_marginal(ll_mat)
  
  waic_res <- suppressWarnings( waic(ll_mat) )
  WAIC_val <- waic_res$estimates["waic","Estimate"]
  
  #BIC computation
  #Sum of the log-likelihoods for each draw
  ll_sum <- rowSums(ll_mat)
  
  # maximum of the sum: estimate of log p(y|theta_hat), that is, log-lik_max
  lhat <- max(ll_sum)
  
  # BIC_val <- -2 * lhat + k * log(N_train)
  BIC_val <- -2 * lhat + k * log(N_train - output$mean$train_metrics_skip)
  
  # BIC Approx
  # BIC_val  <- DIC_val + k * (log(N_train) - 2) 
  
  ############################################################
  #################### OUT-OF-SAMPLE #########################
  ############################################################

  y_pred <- if (is.null(output$mean$yp_onestep_out)) output$mean$yp_out else output$mean$yp_onestep_out

  y_true <- series[(N_train+1):length(series)]
  e <- y_pred - y_true
  
  # MSE
  MSE_val <- mean((e)^2)
  
  #MAE
  MAE_val <- mean(abs(e))
  MAPE_val <- mean(abs(e / y_true)) * 100
  
  #SMAPE
  SMAPE_val <- mean(2 * abs(e) / (abs(y_true) + abs(y_pred))) * 100
  
  #MASE
  scale <- mean(abs(diff(series[1:N_train])))
  MASE_val <- MAE_val / scale
  
 list(
    DIC   = round(DIC_val, 1),
    WAIC  = round(WAIC_val, 1),
    BIC   = round(BIC_val, 1),
    MSE   = round(MSE_val, 3),
    MAE   = round(MAE_val, 3),
    MAPE  = round(MAPE_val, 1),
    SMAPE = round(SMAPE_val, 1),
    MASE  = round(MASE_val, 3)
  )
}

compute_dic_marginal <- function(ll_mat) {
  dev      <- -2 * ll_mat
  dev_sum  <- rowSums(dev)
  Dbar     <- mean(dev_sum)
  pD       <- 0.5 * var(dev_sum)
  DIC      <- Dbar + pD
  return(DIC)
}

# Appends metrics row to provided metrics_df and returns it
# If no arguments are passed, generates and returns a new empty metrics_df
metrics_df_append <- function(metrics_df = NA, new_mets = NA, name="...") {
  if (anyNA(metrics_df)) {
    return(data.frame(
      Model = character(),
      DIC   = numeric(),  
      WAIC  = numeric(), 
      BIC   = numeric(),
      MSE   = numeric(),  
      MAE   = numeric(), 
      MAPE  = numeric(),
      SMAPE = numeric(),  
      MASE  = numeric(),
      stringsAsFactors = FALSE
    ))
    
    
  } else if (!anyNA(metrics_df) && !all(is.na(new_mets))) {
    
    return(rbind(
      metrics_df,
      data.frame(
        Model = name,
        DIC   = new_mets$DIC,
        WAIC  = new_mets$WAIC,
        BIC   = new_mets$BIC,
        MSE   = new_mets$MSE,
        MAE   = new_mets$MAE,
        MAPE  = new_mets$MAPE,
        SMAPE = new_mets$SMAPE,
        MASE  = new_mets$MASE,
        stringsAsFactors = FALSE
      ))
    )
  }
  
  return(metrics_df)
}

JAGS model strings

At this point, we define all our JAGS model strings so as to be able to reuse them multiple times through out the script. We start with with defining the model for a generic \(\text{AR}(n)\) sequence. The auto-regressive order \(n\) is passed to the model by means of the order variable.

\[ y(t) = \left[ \sum_{i=1}^{n} \alpha_i y(t-i) \right] + \eta(t) \qquad \eta \sim \mathcal{N}(0, \sigma^2) \]

For any order \(n ≥ 1\) the following model will correctly implement an \(\text{AR}(n)\) process, which allows us to automatically test multiple AR processes easily and at once.

# AR(n) string
modelAR_string <- "
model {
  ############################################################
  ######################## LIKELIHOOD ########################
  ############################################################
  
  # 1. Bootstrapping the AR process
  mu[1:order]     <- y[1:order]
  yp_in[1:order]  <- y[1:order]
  
  # 2. Recursive step of the AR process
  for (t in (order+1):N_train) {
    mu[t]     <- inprod(alpha_rev, y[(t-order):(t-1)]) + m0
    y[t]      ~  dnorm(mu[t], tau)
    yp_in[t]  ~  dnorm(mu[t], tau)
  }
  
  for (t in (train_metrics_skip+1):N_train) {
    loglik[t-train_metrics_skip] <- logdensity.norm(y[t], mu[t], tau)
  }
  
  ############################################################
  ######################## PREDICTION ########################
  ############################################################
  
  # Init helper array z
  z[1:order] <- y[(N_train-order+1):N_train]
  
  for (t in (order+1):(order+N_test)) {
    mu_z[t] <- inprod(alpha_rev, z[(t-order):(t-1)]) + m0
    z[t]    ~  dnorm(mu_z[t], tau)
  }
  
  yp_out <- z[(order+1):(order+N_test)]
  
  ############################################################
  ############ ONE STEP AHEAD PREDICTION #####################
  ############################################################
  
  for (t in 1:N_test) {
    mu_onestep[t]     <- inprod(alpha_rev, y[(N_train+t-order):(N_train+t-1)]) + m0
    yp_onestep_out[t]  ~ dnorm(mu_onestep[t], tau)
  }
  
  ############################################################
  ########################## PRIOR ###########################
  ############################################################
  
  m0  ~ dnorm(0.0, 1.0E-4)
  tau ~ dgamma(0.1, 10)
  for (i in 1:order) {
    alpha[i] ~ dunif(-1.5, 1.5)
    alpha_rev[order - i + 1] <- alpha[i]
  }
  
  # Defining converstion from precision to variance here
  sigma2 <- 1 / tau
}"

We now define a generic \(\text{MA}(n)\) process. The memory order \(n\) is passed to the model by means of the order variable.

\[ y(t) = \left[ \sum_{i=1}^n \beta_i \eta(t-i) \right] + \eta(t) \qquad \eta \sim \mathcal{N}(0, \sigma^2) \]

For any order \(n ≥ 1\) the following model will correctly implement a \(\text{MA}(n)\) process, which allows us to automatically test multiple MA processes easily and at once.

modelMA_string <- "
model {

  ############################################################
  ######################## LIKELIHOOD ########################
  ############################################################

  e[1:order] = rep(0, order)
  yp_in[1:order]  <- y[1:order]
  
  # 2. Recursive step of the MA process
  for (t in (order+1):N_train) {
    mu[t]    <- m0 + inprod(beta_rev, e[(t-order):(t-1)])
    y[t]      ~ dnorm(mu[t], tau)
    yp_in[t]  ~ dnorm(mu[t], tau)
    e[t]     <- y[t] - mu[t]
  }

  for (t in (train_metrics_skip+1):N_train) {
    loglik[t-train_metrics_skip] <- logdensity.norm(y[t], mu[t], tau)
  }
  
  ############################################################
  ######################## PREDICTION ########################
  ############################################################

  for (t in 1:N_test) {
    e[N_train+t] ~ dnorm(0, tau)
    mu_out[t]   <- m0 + inprod(beta_rev, e[(N_train+t-order):(N_train+t-1)])
    yp_out[t]    ~ dnorm(mu_out[t], tau)
  }
  
  ############################################################
  ############ ONE STEP AHEAD PREDICTION #####################
  ############################################################
  
  e_onestep[1:order]   <- e[(N_train-order+1):N_train]
  for (t in 1:N_test) {
    mu_onestep[t]      <- m0 + inprod(beta_rev, e_onestep[t:(t+order-1)])
    yp_onestep_out[t]   ~ dnorm(mu_onestep[t], tau)
    e_onestep[t+order] <- y[N_train+t] - mu_onestep[t]
  }
  
  ############################################################
  ########################## PRIOR ###########################
  ############################################################
  
  m0    ~ dnorm(0.0, 1.0E-4)
  tau   ~ dgamma(0.01, 0.01)

  for (i in 1:order){
    beta[i]      ~ dnorm(0, 1/10)
    beta_rev[i] <- beta[order-i+1]
  }
  
  # Defining converstion from precision to variance here
  sigma2 <- 1 / tau
}"

Finally, the following model will merge the two previous models defining an \(\text{ARMA}(p, q)\) model.

\[ y(t) = \left[ \sum_{i=1}^p \alpha_i y(t-i) \right] + \left[ \sum_{i=1}^q \beta_i \eta(t-i) \right] + \eta(t) \qquad \eta \sim \mathcal{N}(0, \sigma^2). \]

The parameters order_ar and order_ma allow us to specify the two orders when simulating the model so as to be able to test all possible combinations of auto-regressive and moving-average processes.

modelARMA_string <- "
model {
  
  ############################################################
  ######################## LIKELIHOOD ########################
  ############################################################

  e[1:order]     <- rep(0, order)
  mu[1:order]    <- y[1:order]
  yp_in[1:order] <- y[1:order]
  
  # 2. Recursive step of the MA process
  for (t in (order+1):N_train) {
    mu_ar[t] <- inprod(alpha_rev, y[(t-order):(t-1)])
    mu_ma[t] <- inprod(beta_rev, e[(t-order):(t-1)])
    mu[t]    <- m0 + mu_ar[t] + mu_ma[t]
    y[t]      ~ dnorm(mu[t], tau)
    yp_in[t]  ~ dnorm(mu[t], tau)
    e[t]     <- y[t] - mu[t]
  }

  for (t in (train_metrics_skip+1):N_train) {
    loglik[t-train_metrics_skip] <- logdensity.norm(y[t], mu[t], tau)
  }
  
  ############################################################
  ######################## PREDICTION ########################
  ############################################################

  z[1:order]   <- y[(N_train-order+1):N_train]
  e_z[1:order] <- e[(N_train-order+1):N_train]

  for (t in (order+1):(order+N_test)) {
    e_z[t]        ~ dnorm(0, tau)
    mu_ar_out[t] <- inprod(alpha_rev, z[(t-order):(t-1)])
    mu_ma_out[t] <- inprod(beta_rev, e_z[(t-order):(t-1)])
    mu_z[t]      <- m0 + mu_ar_out[t] + mu_ma_out[t]
    z[t]          ~ dnorm(mu_z[t], tau)
  }

  yp_out <- z[(order+1):(order+N_test)]
  
  ############################################################
  ############ ONE STEP AHEAD PREDICTION #####################
  ############################################################
  
  e_onestep[1:order]   <- e[(N_train-order+1):N_train]
  for (t in 1:N_test) {
    mu_ar_onestep[t] <- inprod(alpha_rev, y[(N_train+t-order):(N_train+t-1)])
    mu_ma_onestep[t]   <- inprod(beta_rev, e_onestep[t:(t+order-1)])
    mu_onestep[t]      <- m0 + mu_ar_onestep[t] + mu_ma_onestep[t]
    yp_onestep_out[t]   ~ dnorm(mu_onestep[t], tau)
    e_onestep[t+order] <- y[N_train+t] - mu_onestep[t]
  }
  
  ############################################################
  ########################## PRIOR ###########################
  ############################################################
  
  m0    ~ dnorm(0.0, 1.0E-4)
  tau   ~ dgamma(0.01, 0.01)
  
  for (i in 1:order_ar) {
    alpha[i] ~ dunif(-1.5, 1.5)
  }
  for (i in (order_ar+1):order) {
    alpha[i]      <- 0
  }
  for (i in 1:order) {
    alpha_rev[i] <- alpha[order-i+1]
  }

  for (i in 1:order_ma) {
    beta[i] ~ dnorm(0, 1/10)
  }
  for (i in (order_ma+1):order) {
    beta[i]      <- 0
  }
  for (i in 1:order) {
    beta_rev[i] <- beta[order-i+1]
  }
  
  # Defining converstion from precision to variance here
  sigma2 <- 1 / tau
  
}"

Modelling

At this point we can start modelling our two time series using the stochastic processes for which we defined the previous model strings. We will start with the GDP time-series and then focus on the CPI time-series. Finally, we will also attempt to model a correlated process using a \(\text{VAR}(n)\) model, where both time-series will be modelled together at once.

GDP time series

# Metrics dataframe setup
metrics_df_GDP <- metrics_df_append()

We will, at first, start off with some \(\text{AR}(n)\) models, with increasing values of \(n\). We will use the obtained results ot decide whether or not the selected model may apply to the problem at hand, and we will also evaluate which regressive order is best.

to_save <- c(
  "alpha",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our AR(n) from n=1 to max_order
max_order <- if (intense) 9 else 2
for (order in 1:max_order) {
  data_list <- list(
    "y"=series_GDP,
    "order"=order,
    "N_train"=N_train,
    "N_test"=N_test,
    "train_metrics_skip"=50
  )
  
  # Perform simulations and visualise output
  output <- fit_JAGS(
    modelAR_string,
    data_list,
    to_save
  )
  visualise(x_axis, series_GDP, output, title=glue("GDP series, AR({order})"))
  
  # Update metrics
  mets <- compute_metrics(output, series_GDP, N_train, k=order+2)
  metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets, glue("AR({order})"))
}

We will now also fit some \(\text{MA}(n)\) processes (for now \(1 \leq n \leq 5\)) to evaluate whether anyone might be better.

to_save <- c(
  "beta",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our MA(n) from n=1 to max_order
max_order <- if (intense) 5 else 2
for (order in 1:max_order) {
  data_list <- list(
    "y"=series_GDP,
    "order"=order,
    "N_train"=N_train,
    "N_test"=N_test,
    "train_metrics_skip"=50
  )
  
  # Perform simulations and visualise output
  output <- fit_JAGS(
    modelMA_string,
    data_list,
    to_save
  )
  visualise(x_axis, series_GDP, output, title=glue("GDP series, MA({order})"))
  
  # Updating metrics
  mets <- compute_metrics(output, series_GDP, N_train, k=order+2)
  metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets, glue("MA({order})"))
}

Finally, we combine the auto-regressive and moving-average proccesses and try to model our series using an \(\text{ARMA}(p, q)\) process. For now, we will attempt modelling for all combinations of \(1 \leq p \leq 4\) and \(1 \leq q \leq 4\).

to_save <- c(
  "alpha",
  "beta",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our AR(n) from n=1 to max_order
max_order_ar <- if (intense) 4 else 2
max_order_ma <- if (intense) 4 else 2
for (order_ma in 1:max_order_ma) {
  for (order_ar in 1:max_order_ar) {
    data_list <- list(
      "y"=series_GDP,
      "N_train"=N_train,
      "N_test"=N_test,
      "order"=max(order_ar, order_ma),
      "order_ar"=order_ar,
      "order_ma"=order_ma,
      "train_metrics_skip"=50
    )
    
    # Perform simulations and visualise output
    output <- fit_JAGS(
      modelARMA_string,
      data_list,
      to_save
    )
    visualise(x_axis, series_GDP, output, title=glue("GDP series, ARMA({order_ar}, {order_ma})"))
    
    # Updating metrics
    mets <- compute_metrics(output, series_GDP, N_train, k=order_ar+order_ma+2)
    metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets, glue("ARMA({order_ar}, {order_ma})"))
  }
}

CPI time series

# Metrics dataframe setup
metrics_df_CPI <- metrics_df_append()

We continue our modelling moving to the CPI time-series. We will proceede as for the GDP, trying our a few \(\text{AR}(n)\) models, then moving to some \(\text{AR}(n)\) and finally, putting both together, some \(\text{ARMA}(p, q)\) models.

to_save <- c(
  "alpha",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our AR(n) from n=1 to max_order
max_order <- if (intense) 9 else 2
for (order in 1:max_order) {
  data_list <- list(
    "y"=series_CPI,
    "order"=order,
    "N_train"=N_train,
    "N_test"=N_test,
    "train_metrics_skip"=50
  )
  
  # Perform simulations and visualise output
  output <- fit_JAGS(
    modelAR_string,
    data_list,
    to_save
  )
  visualise(x_axis, series_CPI, output, title=glue("CPI series, AR({order})"))
  
  if (order == 9) {
    output_ar9_cpi = output
  }
  
  # Update metrics
  mets <- compute_metrics(output, series_CPI, N_train, k=order+2)
  metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets, glue("AR(",order,")"))
}

We will now also fit some \({MA}(n)\) processes (for now \(1 \leq n \leq 5\)) to evaluate whether anyone might be better.

to_save <- c(
  "beta",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our MA(n) from n=1 to max_order
max_order <- if (intense) 5 else 2
for (order in 1:max_order) {
  data_list <- list(
    "y"=series_CPI,
    "order"=order,
    "N_train"=N_train,
    "N_test"=N_test,
    "train_metrics_skip"=50
  )
  
  # Perform simulations and visualise output
  output <- fit_JAGS(
    modelMA_string,
    data_list,
    to_save
  )
  visualise(x_axis, series_CPI, output, title=glue("CPI series, MA({order})"))
  
  # Updating metrics
  mets <- compute_metrics(output, series_CPI, N_train, k=order+2)
  metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets, glue("MA(",order,")"))
}

Finally, we combine the auto-regressive and moving-average proccesses and try to model our series using an \(\text{ARMA}(p, q)\) process. For now, we will attempt modelling for all combinations of \(1 \leq p \leq 4\) and \(1 \leq q \leq 4\).

to_save <- c(
  "alpha",
  "beta",
  "m0",
  "sigma2",
  "yp_in",
  "yp_out",
  "yp_onestep_out",
  "train_metrics_skip",
  "loglik"
)

# Trying our AR(n) from n=1 to max_order
max_order_ar <- if (intense) 4 else 2
max_order_ma <- if (intense) 4 else 2
for (order_ma in 1:max_order_ma) {
  for (order_ar in 1:max_order_ar) {
    data_list <- list(
      "y"=series_CPI,
      "N_train"=N_train,
      "N_test"=N_test,
      "order"=max(order_ar, order_ma),
      "order_ar"=order_ar,
      "order_ma"=order_ma,
      "train_metrics_skip"=50
    )
    
    # Perform simulations and visualise output
    output <- fit_JAGS(
      modelARMA_string,
      data_list,
      to_save
    )
    visualise(x_axis, series_CPI, output, title=glue("CPI series, ARMA({order_ar}, {order_ma})"))
    
    # Updating metrics
    mets <- compute_metrics(output, series_CPI, N_train, k=order_ar+order_ma+2)
    metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets, glue("ARMA({order_ar}, {order_ma})"))
  }
}

VAR modelling

Finally, we will try to fit generic order \(\text{VAR}(n)\) models, which have the following form:

\[ \vec{y}(t) = \vec{\mu}_0 + \left[ \sum_{i=1}^n A_i \vec{y}(t-i) \right] + \vec{\eta}(t) \qquad \forall i \; A_i \in \mathbb{M}(2 \times 2, \mathbb{R}) \qquad \vec{\eta} \sim \mathcal{N}(\vec{0}, \Sigma) \]

Where \(\Sigma\) is the covariance matrix of our white noise (now an \(\mathbb{R}^2\) vector) and \(\vec{y}(t) \in \mathbb{R}^2\) now contains both the GDP and CPI time-series, which will be modelled together in a correlated fashion. By passing the parameter order during simulation we can specify the auto-regressive order \(n\) of the multivariate process.

The idea behind this model is to attempt to model not only the two time-series independently, but also in a correlated way, as we work under the assumption that the two may influence each other.

# Prior fragment for Wishart distributions (positive definite covariance)
prior_wishart <- "
  # Wishart prior for a full covariance matrix
  R[1,1] <- 1
  R[2,2] <- 1
  R[1,2] <- 0
  R[2,1] <- 0
  tau ~ dwish(R, 3) # df=3 is the minimum for a 2x2 matrix
"

# Prior fragment for independent Gamma distributions (diagonal covariance)
prior_diagonal <- "
  # Independent Gamma priors for diagonal elements of precision matrix
  tau[1,1] ~ dgamma(0.01, 0.01)
  tau[2,2] ~ dgamma(0.01, 0.01)
  tau[1,2] <- 0 
  tau[2,1] <- 0
"

modelVAR_string = "
model {

  ############################################################
  ######################## LIKELIHOOD ########################
  ############################################################
  
  mu[1:2,1:order] = y[,1:order]
  yp_in[1:2,1:order] = y[,1:order]
  
  for (t in (order+1):N_train) {
    mu_tmp[1:2, t, 1] <- m0
    for (j in 1:order) {
      mu_tmp[1:2,t,j+1] <- mu_tmp[1:2, t, j] + A[,,j] %*% y[,(t-j)]
    }
    mu[1:2,t]   <- mu_tmp[1:2,t,order+1]
    y[1:2,t]     ~ dmnorm(mu[1:2,t], tau)
    yp_in[1:2,t] ~ dmnorm(mu[1:2,t], tau)
  }

  for (t in (train_metrics_skip+1):N_train) {
    loglik_gdp[t-train_metrics_skip]  <- logdensity.norm(y[1,t], mu[1,t], 1 / sigma2[1,1])
    loglik_infl[t-train_metrics_skip] <- logdensity.norm(y[2,t], mu[2,t], 1 / sigma2[2,2])
  }

  ############################################################
  ######################## PREDICTION ########################
  ############################################################
  
  z[1:2,1:order] = y[,(N_train-order+1):N_train]
  
  for (t in (order+1):(order+N_test)) {
    mu_z[1:2, t, 1] <- m0
    for (j in 1:order) {
      mu_z[1:2,t,j+1] <- mu_z[1:2, t, j] + A[,,j] %*% z[,(t-j)]
    }
    z[1:2,t] ~ dmnorm(mu_z[1:2,t,order+1], tau)
  }
  
  yp_out[1:2,1:N_test] = z[1:2,(order+1):(order+N_test)]
  
  ############################################################
  ############ ONE STEP AHEAD PREDICTION #####################
  ############################################################
  
  for (t in 1:N_test) {
    mu_onestep_tmp[1:2, t, 1] <- m0
    for (j in 1:order) {
      mu_onestep_tmp[1:2,t,j+1] <- mu_onestep_tmp[1:2, t, j] + A[,,j] %*% y[,(N_train+t-j)]
    }
    mu_onestep[1:2,t]     <- mu_onestep_tmp[1:2,t,order+1]
    yp_onestep_out[1:2,t]  ~ dmnorm(mu_onestep[1:2,t], tau)
  }
  
  ############################################################
  ########################## PRIOR ###########################
  ############################################################

  for (i in 1:2) {
    for (j in 1:2) {
      for (o in 1:order) {
        A[i, j, o] ~ dunif(-1.2, 1.2)
        # A[i, j, o] ~ dnorm(0, 10)
      }
    }
  }
  
  # Dynamically inserted prior for tau
  
  {{prior_for_tau}}
  
  # __if UNCORRELATED_ERRORS__
  # 
  # tau[1,1] ~ dgamma(0.01, 0.01)
  # tau[2,2] ~ dgamma(0.01, 0.01)
  # tau[1,2] <- 0
  # tau[2,1] <- 0
  
  # __else__
  # 
  # tau     ~ dwish(R, 3) # dwish(R=ScaleMatrix... Identity, df=degrees of freedom)
  
  m_gdp   ~ dnorm(0.0, 1.0E-4)
  m_infl  ~ dnorm(0.0, 1.0E-4)
  m0     <- c(m_gdp, m_infl)
  
  sigma2 <- inverse(tau)
  
  ############################################################
  ######################## EXTRACTION ########################
  ############################################################

  # Expose them for output in R
  loglik_gdp_out   <- loglik_gdp
  loglik_infl_out  <- loglik_infl

  yp_in_gdp   <- yp_in[1,]
  yp_in_infl  <- yp_in[2,]
  yp_out_gdp  <- yp_out[1,]
  yp_out_infl <- yp_out[2,]
  yp_onestep_out_gdp  <- yp_onestep_out[1,]
  yp_onestep_out_infl <- yp_onestep_out[2,]
}"

# Setting up parameters
y_var <- matrix(c(series_GDP, series_CPI), nrow=2, byrow=TRUE)
order <- 1
use_wishart <- TRUE

if (use_wishart) {
  prior_for_tau = prior_wishart
} else {
  prior_for_tau = prior_diagonal
}
final_modelVAR_string <- glue(modelVAR_string, .open="{{", .close="}}")
# Uncomment for debugging inspection
#print(glue(modelVAR_string, .open="{{", .close="}}"))

# Defining JAGS parameters
to_save <- c(
  "A",
  "m_gdp",
  "m_infl",
  "sigma2",
  "yp_in_gdp",
  "yp_in_infl",
  "yp_out_gdp",
  "yp_out_infl",
  "yp_onestep_out_gdp",
  "yp_onestep_out_infl",
  "train_metrics_skip",
  "loglik_gdp",   # marginal log-likelihood for GDP
  "loglik_infl"   # marginal log-likelihood for CPI
)
data_list <- list(
  "y"=y_var,
  "order"=order,
  "N_train" = N_train,
  "N_test" = N_test,
  "train_metrics_skip"=50
)

# Fitting VAR model and visualising it
jags_output_VAR = fit_JAGS(
  final_modelVAR_string,
  data_list,
  to_save,
  n_iter=6000,
  n_adapt=1000,
  n_chain=1,
  n_burnin=2000
)
visualise_VAR(
  x_axis,
  series_GDP,
  series_CPI,
  jags_output_VAR,
  title=glue("VAR({order})"),
  keep=5
)

# Wrapper for GDP-only metrics
output_gdp <- list(
  mean      = list(
    yp_out = jags_output_VAR$mean$yp_out_gdp,
    yp_onestep_out = jags_output_VAR$mean$yp_onestep_out_gdp,
    train_metrics_skip=jags_output_VAR$mean$train_metrics_skip
  ),
  sims.list = list(loglik  = jags_output_VAR$sims.list$loglik_gdp)
)

# Wrapper for CPI-only metrics
output_infl <- list(
  mean      = list(
    yp_out = jags_output_VAR$mean$yp_out_infl,
    yp_onestep_out = jags_output_VAR$mean$yp_onestep_out_infl,
    train_metrics_skip=jags_output_VAR$mean$train_metrics_skip
  ),
  sims.list = list(loglik  = jags_output_VAR$sims.list$loglik_infl)
)

# Compute univariate metrics for GDP and CPI as well
# RMK: just for visualisation purposes: can't compare marginalised metrics with univariate metrics
mets_VAR_gdp  <- compute_metrics(output_gdp,  series_GDP, N_train, k = 4*order + 2 + 3)
mets_VAR_infl <- compute_metrics(output_infl, series_CPI, N_train, k = 4*order + 2 + 3)

metrics_df_GDP <- metrics_df_append(metrics_df_GDP, mets_VAR_gdp, glue("VAR_univar_GDP({order})"))
metrics_df_CPI <- metrics_df_append(metrics_df_CPI, mets_VAR_infl, glue("VAR_univar_GDP({order})"))

VAR: Correlated vs Uncorrelated errors

We will now try to fit the \(\text{VAR}(1)\) model using two different priors for the precision matrix of \(\tau\) of the 2 errors.

Prior 1: diagonal precision matrix

We will force the errors to be uncorrelated. We will do this by setting the priors for the off-diagonal elements to be 0, while leaving a gamma prior to the diagonal elements, akin to the priors used of the univariate AR models.

Prior 2: positive definite precision matrix

We will set a Wishart prior to the precision matrix of \(\tau\). This will allow the precision matrix to have off-diagonal elements different from zero, but will still ensure that it will be positive definite.

# Uncorrelated errors
prior_for_tau         <- prior_diagonal
modelVAR_string_gamma <- glue(modelVAR_string, .open="{{", .close="}}")
jags_output_VAR_gamma <- fit_JAGS(
  modelVAR_string_gamma,
  data_list,
  to_save,
  n_iter=6000,
  n_adapt=1000,
  n_chain=1,
  n_burnin=2000
)

# Possibly correlated errors (Wishart prior)
prior_for_tau           <- prior_wishart
modelVAR_string_wishart <- glue(modelVAR_string, .open="{{", .close="}}")
jags_output_VAR_wishart <- fit_JAGS(
  modelVAR_string_wishart,
  data_list,
  to_save,
  n_iter=6000,
  n_adapt=1000,
  n_chain=1,
  n_burnin=2000
)

# Logging results of the two models
print(glue("Posterior means using diagonal prior"))

## Posterior means using diagonal prior

print(glue("\r\nA"))

## 
## A

print(jags_output_VAR_gamma$mean$A[,,1])

##            [,1]        [,2]
## [1,] 0.90607191 -0.03350507
## [2,] 0.08042365  0.88547526

print(glue("\r\nsigma2"))

## 
## sigma2

print(jags_output_VAR_gamma$mean$sigma2)

##          [,1]      [,2]
## [1,] 2.607159 0.0000000
## [2,] 0.000000 0.9587027

print(glue("\r\n\r\nPosterior means using Wishart prior"))

## 
## 
## Posterior means using Wishart prior

print(glue("\r\nA"))

## 
## A

print(jags_output_VAR_wishart$mean$A[,,1])

##            [,1]        [,2]
## [1,] 0.90977027 -0.03879845
## [2,] 0.07861106  0.88670588

print(glue("\r\nsigma2"))

## 
## sigma2

print(jags_output_VAR_wishart$mean$sigma2)

##           [,1]      [,2]
## [1,] 2.5826247 0.5496647
## [2,] 0.5496647 0.9546752

From these results we see that the noises of the 2 time series are somewhat correlated. We can compute the Pearson Correlation Coefficient to further inspect this hypotheses:

\[ \rho(err_{gdp},err_{cpi}) \;=\;\frac{\text{Cov}(err_{gdp},err_{cpi})}{\sqrt{\text{Var}(err_{gdp})\,\times\,\text{Var}(err_{cpi})}} \;=\; \frac{0.5529124}{\sqrt{2.5870927 \times 0.9520914}} \;\approx\;0.3523 \]

This indicates a moderate positive linear relationship between the two variables. It suggests that as one variable increases, the other tends to increase as well, but not perfectly.

On the other hand, the effects of the 2 different priors on the matrix A is very moderate.

This suggest that the two variables are moderately related through their contemporaneous shocks, but this relation is not dominant enough to significantly affect the lagged relationship between the two variables.

Metrics

At this point, we display all so far computer metrics for all models and both time-series, so as to be able to perform model selection objectively.

# GDP table
kable(
  metrics_df_GDP,
  caption = "Univariate model metrics — GDP series",
  digits  = c(1,1,1,1,3,3,1,1,3),
  format  = "html",
  table.attr = 'cellpadding="10" cellspacing="0"'
) %>%
  kable_styling(full_width = F, position = "center") %>%
  row_spec(0, background = "#000000", color = "#FFFFFF") %>%  # intestazione
  row_spec(1:nrow(metrics_df_GDP), background = "#000000", color = "#FFFFFF")  # righe dati

Univariate model metrics — GDP series

Model	DIC	WAIC	BIC	MSE	MAE	MAPE	SMAPE	MASE
AR(1)	744.7	705.6	694.4	7.971	1.378	75.0	26.5	1.195
AR(2)	690.0	667.2	666.5	12.174	1.828	144.0	31.9	1.585
AR(3)	690.5	666.3	666.8	11.280	1.731	134.9	30.5	1.501
AR(4)	682.5	663.8	673.2	11.885	1.756	135.7	30.0	1.523
AR(5)	644.2	630.5	647.0	8.451	1.543	124.5	27.0	1.338
AR(6)	647.2	632.9	653.0	8.151	1.522	121.9	26.9	1.320
AR(7)	654.9	636.7	658.0	8.193	1.532	121.0	25.9	1.329
AR(8)	660.6	637.8	660.1	8.373	1.590	125.6	28.3	1.379
AR(9)	645.4	626.7	655.6	7.683	1.494	124.1	26.6	1.296
MA(1)	846.8	835.5	838.8	9.079	2.077	114.2	42.6	1.801
MA(2)	792.0	780.5	786.8	15.539	2.488	130.8	45.0	2.157
MA(3)	659.1	638.7	643.1	5.097	1.316	100.7	27.8	1.142
MA(4)	620.8	604.1	613.0	5.912	1.344	140.0	27.6	1.166
MA(5)	652.9	620.5	627.9	7.048	1.388	133.0	25.5	1.203
ARMA(1, 1)	697.9	673.7	672.0	11.378	1.880	136.7	34.7	1.630
ARMA(2, 1)	684.6	664.0	669.9	11.637	1.789	138.8	31.4	1.551
ARMA(3, 1)	675.6	654.3	662.1	11.163	1.799	165.1	32.8	1.560
ARMA(4, 1)	690.3	661.6	667.4	9.795	1.701	141.0	31.9	1.475
ARMA(1, 2)	683.2	671.7	661.2	6.573	1.457	115.2	28.5	1.263
ARMA(2, 2)	646.4	631.6	639.4	7.314	1.559	136.2	30.5	1.352
ARMA(3, 2)	632.2	613.6	626.9	7.601	1.567	132.5	30.3	1.359
ARMA(4, 2)	649.3	624.2	638.6	7.902	1.560	129.1	30.4	1.353
ARMA(1, 3)	592.2	585.0	598.6	6.245	1.295	142.4	24.9	1.123
ARMA(2, 3)	598.2	571.4	582.2	5.323	1.099	130.7	22.1	0.953
ARMA(3, 3)	594.0	579.7	598.5	5.837	1.204	137.2	23.5	1.044
ARMA(4, 3)	611.2	590.1	610.5	5.958	1.251	136.6	24.2	1.085
ARMA(1, 4)	601.2	582.0	595.5	6.111	1.260	138.7	24.2	1.093
ARMA(2, 4)	614.7	583.4	593.7	6.259	1.285	142.0	24.8	1.115
ARMA(3, 4)	608.7	587.2	606.7	6.126	1.267	138.2	24.4	1.098
ARMA(4, 4)	611.3	589.3	615.6	5.957	1.253	136.9	24.3	1.086
VAR_univar_GDP(1)	740.9	704.1	725.0	8.013	1.381	76.1	26.9	1.198

# CPI table
kable(
  metrics_df_CPI,
  caption = "Univariate model metrics — CPI series",
  digits  = c(1,1,1,1,3,3,1,1,3),
  format  = "html",
  table.attr = 'cellpadding="10" cellspacing="0"'
) %>%
  kable_styling(full_width = F, position = "center") %>%
  row_spec(0, background = "#000000", color = "#FFFFFF") %>%
  row_spec(1:nrow(metrics_df_CPI), background = "#000000", color = "#FFFFFF")

Univariate model metrics — CPI series

Model	DIC	WAIC	BIC	MSE	MAE	MAPE	SMAPE	MASE
AR(1)	569.2	552.4	550.4	0.689	0.607	202.3	33.5	0.899
AR(2)	562.1	548.3	552.0	0.580	0.592	217.8	35.0	0.877
AR(3)	557.6	546.6	554.1	0.559	0.578	200.9	37.8	0.857
AR(4)	542.0	532.6	546.4	0.566	0.571	169.0	36.4	0.847
AR(5)	509.7	502.2	520.3	0.508	0.575	152.1	36.4	0.853
AR(6)	506.9	498.2	518.2	0.511	0.580	183.8	36.5	0.860
AR(7)	510.5	499.8	520.9	0.508	0.577	191.6	36.3	0.855
AR(8)	511.2	500.4	526.3	0.515	0.581	190.2	36.0	0.861
AR(9)	482.6	468.2	499.0	0.447	0.528	263.3	34.3	0.783
MA(1)	800.3	792.5	797.2	1.887	1.181	913.3	52.7	1.752
MA(2)	737.9	725.9	727.1	1.321	0.972	651.8	47.1	1.441
MA(3)	571.7	560.5	571.4	0.589	0.594	441.3	34.0	0.880
MA(4)	526.1	507.8	515.0	0.573	0.579	423.0	35.5	0.859
MA(5)	551.3	526.7	536.1	0.467	0.531	336.0	34.9	0.787
ARMA(1, 1)	556.3	546.8	551.4	0.606	0.605	217.9	35.9	0.898
ARMA(2, 1)	546.5	527.7	527.5	0.619	0.632	276.1	41.3	0.937
ARMA(3, 1)	546.7	540.8	556.6	0.561	0.577	201.1	37.3	0.855
ARMA(4, 1)	511.5	499.6	509.1	0.537	0.609	268.6	40.4	0.903
ARMA(1, 2)	541.8	532.4	534.1	0.529	0.524	143.7	34.5	0.777
ARMA(2, 2)	543.3	526.2	531.4	0.638	0.639	304.3	41.6	0.948
ARMA(3, 2)	501.7	494.5	509.2	0.473	0.508	141.1	35.0	0.753
ARMA(4, 2)	504.5	497.5	514.1	0.566	0.605	291.2	40.3	0.898
ARMA(1, 3)	477.8	468.0	479.2	0.424	0.534	304.3	35.0	0.792
ARMA(2, 3)	467.7	456.5	471.7	0.401	0.520	259.5	34.8	0.770
ARMA(3, 3)	436.9	421.3	435.7	0.378	0.496	230.3	35.0	0.735
ARMA(4, 3)	467.2	447.7	464.2	0.436	0.525	266.5	36.0	0.779
ARMA(1, 4)	464.1	449.1	462.4	0.405	0.493	222.1	33.5	0.731
ARMA(2, 4)	478.8	449.0	458.2	0.426	0.484	221.6	33.2	0.717
ARMA(3, 4)	476.6	452.2	462.1	0.410	0.489	219.6	33.5	0.725
ARMA(4, 4)	457.2	440.5	464.1	0.449	0.525	214.4	36.4	0.779
VAR_univar_GDP(1)	550.5	537.6	567.1	0.721	0.636	200.4	39.3	0.944

NOTICE: The MAPE metric exceeds 100% because the series (quarterly or monthly GDP growth rates) has true values (denominators) that are very small or close to zero. This results in an amplified percentage.

Trace plots and posterior densities

We will now display some trace plots for some parameters of some of the fitted models and relative posterior density plots to make sure the MCMC samples are not getting stuck in some local regions, and are behaving well overall.

AR(9) trace and density plots (GDP)

n_chain = if (intense) 3 else 1

to_save <- c(
  "alpha",
  "m0",
  "sigma2"
)

data_list <- list(
  "y"=series_GDP,
  "order"=9,
  "N_train"=N_train,
  "N_test"=N_test,
  "train_metrics_skip"=50
)
  
output <- fit_JAGS(
  modelAR_string,
  data_list,
  to_save,
  n_thin=100,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=100000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 251
##    Unobserved stochastic nodes: 352
##    Total graph size: 1901
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 99000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

color_scheme_set("teal")
mcmc_trace(output$samples[,1:11])

mcmc_dens(output$samples[,1:11])

mcmc_intervals(output$samples[,1:11])

MA(3) trace and density plots (GDP)

n_chain = if (intense) 3 else 1

to_save <- c(
  "beta",
  "m0",
  "sigma2"
)

data_list <- list(
  "y"=series_GDP,
  "order"=3,
  "N_train"=N_train,
  "N_test"=N_test,
  "train_metrics_skip"=50
)
  
output <- fit_JAGS(
  modelMA_string,
  data_list,
  to_save,
  n_thin=100,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=100000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 257
##    Unobserved stochastic nodes: 397
##    Total graph size: 2265
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 99000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

mcmc_trace(output$samples[,1:5])

mcmc_dens(output$samples[,1:5])

mcmc_intervals(output$samples[,1:5])

ARMA order selection and traces (GDP)

## --- 1. Support Function: matrix [p × q] ---------------------------
make_matrix <- function(df, metric, P, Q){
  M <- matrix(NA_real_, nrow = P, ncol = Q)   # row = p,  column = q
  for(i in seq_len(nrow(df))){
    nums <- as.integer(regmatches(df$Model[i],
                   gregexpr("\\d+", df$Model[i]))[[1]])
    if(length(nums) == 2){
      p <- nums[1];  q <- nums[2]
      M[p, q] <- df[i, metric]
    }
  }
  M
}
## --- 2. function for plotting axis  -------------------------
plot_heat <- function(mat, title){
  image(x = 1:ncol(mat),            
        y = 1:nrow(mat),            
        z = mat,                  
        col = colorRampPalette(
                c("navy","deepskyblue","yellow","orange","red"))(100),
        axes = FALSE, xlab = "", ylab = "", main = title, useRaster = TRUE)
  contour(1:ncol(mat), 1:nrow(mat), mat,
          add = TRUE, lwd = 0.7, drawlabels = FALSE)
  axis(1, at = 1:ncol(mat), labels = 1:ncol(mat), cex.axis = 0.75) # q
  axis(2, at = 1:nrow(mat), labels = 1:nrow(mat), las = 1,
       cex.axis = 0.75)                                           # p
  box()
}

## --- 3. parameters & palette --------------------------------------------
max_p <- max_order_ar   # max order AR
max_q <- max_order_ma   # max order MA
metrics <- c("DIC","WAIC","BIC","MSE",
             "MAE","MAPE","SMAPE","MASE")

## --- 4. heat-map 2×4 for axis p↔q --------------------------------------
op <- par(mfrow = c(2,4), mar = c(2.2,2.2,2.4,0.6))

for(met in metrics){
  mat <- make_matrix(metrics_df_GDP, met, max_p, max_q)
  plot_heat(mat, paste("ARMA –", met))
}

par(op)

n_chain = if (intense) 3 else 1

to_save <- c(
  "alpha",
  "beta",
  "m0",
  "sigma2"
)

order_ar <- 2
order_ma <- 3
data_list <- list(
  "y"=series_GDP,
  "N_train"=N_train,
  "N_test"=N_test,
  "order"=max(order_ar, order_ma),
  "order_ar"=order_ar,
  "order_ma"=order_ma,
  "train_metrics_skip"=50
)

# Perform simulations and visualise output
output <- fit_JAGS(
  modelARMA_string,
  data_list,
  to_save,
  n_thin=50,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=50000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 257
##    Unobserved stochastic nodes: 399
##    Total graph size: 2964
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 49000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

color_scheme_set("teal")
mcmc_trace(output$samples[,1:8])

mcmc_dens(output$samples[,1:8])

mcmc_intervals(output$samples[,1:8])

AR(9) trace and density plots (CPI)

n_chain = if (intense) 3 else 1

to_save <- c(
  "alpha",
  "m0",
  "sigma2"
)

data_list <- list(
  "y"=series_CPI,
  "order"=9,
  "N_train"=N_train,
  "N_test"=N_test,
  "train_metrics_skip"=50
)
  
output <- fit_JAGS(
  modelAR_string,
  data_list,
  to_save,
  n_thin=100,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=100000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 251
##    Unobserved stochastic nodes: 352
##    Total graph size: 1901
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 99000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

color_scheme_set("teal")
mcmc_trace(output$samples[,1:11])

mcmc_dens(output$samples[,1:11])

mcmc_intervals(output$samples[,1:11])

MA(4) trace and density plots (CPI)

n_chain = if (intense) 3 else 1

to_save <- c(
  "beta",
  "m0",
  "sigma2"
)

data_list <- list(
  "y"=series_CPI,
  "order"=4,
  "N_train"=N_train,
  "N_test"=N_test,
  "train_metrics_skip"=50
)
  
output <- fit_JAGS(
  modelMA_string,
  data_list,
  to_save,
  n_thin=100,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=100000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 256
##    Unobserved stochastic nodes: 397
##    Total graph size: 2261
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 99000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

mcmc_trace(output$samples[,1:6])

mcmc_dens(output$samples[,1:6])

mcmc_intervals(output$samples[,1:6])

ARMA order selection and traces (CPI)

## --- 1. Support Function: matrix [p × q] ---------------------------
make_matrix <- function(df, metric, P, Q){
  M <- matrix(NA_real_, nrow = P, ncol = Q)  
  for(i in seq_len(nrow(df))){
    nums <- as.integer(regmatches(df$Model[i],
                   gregexpr("\\d+", df$Model[i]))[[1]])
    if(length(nums) == 2){
      p <- nums[1];  q <- nums[2]
      M[p, q] <- df[i, metric]
    }
  }
  M
}

## --- 2. function for plotting axis  -------------------------
plot_heat <- function(mat, title){
  image(x = 1:ncol(mat),         
        y = 1:nrow(mat),          
        z = mat,                   
        col = colorRampPalette(
                c("navy","deepskyblue","yellow","orange","red"))(100),
        axes = FALSE, xlab = "", ylab = "", main = title, useRaster = TRUE)
  contour(1:ncol(mat), 1:nrow(mat), mat,
          add = TRUE, lwd = 0.7, drawlabels = FALSE)
  axis(1, at = 1:ncol(mat), labels = 1:ncol(mat), cex.axis = 0.75) # q
  axis(2, at = 1:nrow(mat), labels = 1:nrow(mat), las = 1,
       cex.axis = 0.75)                                           # p
  box()
}

## --- 3. parameters & palette --------------------------------------------
max_p <- max_order_ar   # max order AR
max_q <- max_order_ma   # max order MA
metrics <- c("DIC","WAIC","BIC","MSE",
             "MAE","MAPE","SMAPE","MASE")

## --- 4. heat-map 2×4 for axis p↔q --------------------------------------
op <- par(mfrow = c(2,4), mar = c(2.2,2.2,2.4,0.6))

for(met in metrics){
  mat <- make_matrix(metrics_df_CPI, met, max_p, max_q)
  plot_heat(mat, paste("ARMA –", met))
}

par(op)

n_chain = if (intense) 3 else 1

to_save <- c(
  "alpha",
  "beta",
  "m0",
  "sigma2"
)

order_ar <- 3
order_ma <- 3
data_list <- list(
  "y"=series_CPI,
  "N_train"=N_train,
  "N_test"=N_test,
  "order"=max(order_ar, order_ma),
  "order_ar"=order_ar,
  "order_ma"=order_ma,
  "train_metrics_skip"=50
)

# Perform simulations and visualise output
output <- fit_JAGS(
  modelARMA_string,
  data_list,
  to_save,
  n_thin=100,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=100000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 257
##    Unobserved stochastic nodes: 400
##    Total graph size: 2965
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 99000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

color_scheme_set("teal")
mcmc_trace(output$samples[,1:8])

mcmc_dens(output$samples[,1:8])

mcmc_intervals(output$samples[,1:8])

VAR(1) trace and density plots

to_save <- c(
  "A",
  "m0",
  "sigma2"
)

data_list <- list(
  "y"=y_var,
  "order"=1,
  "N_train"=N_train,
  "N_test"=N_test,
  "train_metrics_skip"=50
)
  
output <- fit_JAGS(
  final_modelVAR_string,
  data_list,
  to_save,
  n_thin=3,
  n_adapt=2000,
  n_chain=n_chain,
  n_iter=6000,
  verbose=TRUE
)

## 
## Processing function input....... 
## 
## Done. 
##  
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 259
##    Unobserved stochastic nodes: 356
##    Total graph size: 2735
## 
## Initializing model
## 
## Adaptive phase, 2000 iterations x 3 chains 
## If no progress bar appears JAGS has decided not to adapt 
##  
## 
##  Burn-in phase, 1000 iterations x 3 chains 
##  
## 
## Sampling from joint posterior, 5000 iterations x 3 chains 
##  
## 
## Calculating statistics....... 
## 
## Done.

mcmc_trace(output$samples[,1:10])

mcmc_dens(output$samples[,1:10])

mcmc_intervals(output$samples[,1:10])

Model selection

Based on the calculated metrics, the models that seem to perform best are \(\text{ARMA}(2, 3)\) and \(\text{ARMA}(3, 3)\), both for the GDP and inflation time series.

An honorable mention must be given to the \(\text{AR}(9)\) model for the CPI time series, where the out-of-sample prediction for the next year and a half, given only the in-sample data, is almost perfectly aligned with the actual values, as shown in the following plot.

if (intense) {
  visualise(x_axis, series_CPI, output_ar9_cpi, title=glue("CPI series, AR(9)"))
}

Anomaly detection

In this last section we will examine anomalies and outliers which are present in our two time series. So far, we performed model selection and predictions of future values for the series, however we could clearly see that both the training and test set presented samples which deviated quite heavily from their neighboring values. This can be cause of disturbances during the learning phase and also a reduced accuracy during forecasting.

The objective of this section is to detect said outliers and visualise them so as to clearly show which samples may be problematic during fitting and out-of-sample prediction.

In order to do this, for each timestamp \(t\), we consider the incremental values of our series, each of which is to be treated as an independent random variable:

\[ \Delta_t := y_{t} - y_{t-1} \stackrel{iid}{\sim} \mathcal{N}(\mu, \sigma_t^2) \]

In order to be able to perform anomaly detection, we must ensure that the sequence of incremental values \(\Delta_t\) can be approximated to a random walk. This can be verified by trying to fit an \(\text{AR}(1)\) process to the initial sequence of values \(y_t\) and checking whether the value \(|\alpha|\) is close to 1:

# Data lists for both delta-series
data_list_GDP <- list(
  "y"=series_GDP,
  "train_metrics_skip"=50,
  "order"=1,
  "N_train"=N_tot,
  "N_test"=1          # No need to perform forecasting
)
data_list_CPI <- list(
  "y"=series_CPI,
  "train_metrics_skip"=50,
  "order"=1,
  "N_train"=N_tot,
  "N_test"=1          # No need to perform forecasting
)

# Only need to inspect value of alpha for both series
to_save <- c(
  "alpha"
)

# Fit AR(1) models
output_GDP <- fit_JAGS(
  modelAR_string,
  data_list_GDP,
  to_save
)
output_CPI <- fit_JAGS(
  modelAR_string,
  data_list_CPI,
  to_save
)

# Inspecting values of alpha
print(glue("Alpha value for GDP has posterior mean {output_GDP$mean$alpha} with standard deviation of {output_GDP$sd$alpha}"))

## Alpha value for GDP has posterior mean 0.859521726483674 with standard deviation of 0.030267447634924

print(glue("Alpha value for CPI has posterior mean {output_CPI$mean$alpha} with standard deviation of {output_CPI$sd$alpha}"))

## Alpha value for CPI has posterior mean 0.939175750757215 with standard deviation of 0.0200438333942664

Now that we know that the two series can be approximated reasonably well to random walks (\(\alpha_{GDP} = 0.86\) and \(\alpha_{CPI} = 0.938\)) we proceede with our anomaly detection.

Notice that the random variable \(\Delta_t\) has a fixed mean value \(\mu\), but a time dependant variance \(\sigma_t^2\), which we model as:

\[ \sigma_t^{-2} = \tau_t = \beta_1 + \beta_2 \delta_t \qquad \delta_t \sim Ber(p) \]

We can, therefore, define the following JAGS model-string which will fit the individual values of \(\delta_t\).

anomaly_detection_string <- "
model {

  ############################################################
  ######################## LIKELIHOOD ########################
  ############################################################

  for(t in 1:N_tot) {
    delta[t]  ~  dnorm(mu, tau[t])
    tau[t]    <- beta[1] + beta[2] * d[t]
  }
  
  
  ############################################################
  ########################## PRIOR ###########################
  ############################################################
  
  # Beta coefficients, Bernoulli probability and mean
  beta[1] ~ dexp(0.1)
  beta[2] ~ dexp(0.1)
  p       ~ dbeta(1, 1)
  mu      ~ dnorm(0, 0.01)
  
  # Delta value
  for(t in 1:N_tot){
    d[t] ~ dbern(p)
  }
}"

# Obtaining series of differences
delta_GDP <- series_GDP[2:N_tot] - series_GDP[1:(N_tot-1)]
delta_CPI <- series_CPI[2:N_tot] - series_CPI[1:(N_tot-1)]

# Setting up parameters for JAGS
data_list_GDP <- list(
  "delta"=delta_GDP,
  "N_tot"=N_tot-1
)
data_list_CPI <- list(
  "delta"=delta_CPI,
  "N_tot"=N_tot-1
)
to_save <- c(
  "d",
  "mu",
  "tau"
)

# Perform simulations
anomaly_output_GDP <- fit_JAGS(
  anomaly_detection_string,
  data_list_GDP,
  to_save,
  n_iter=10000
)
anomaly_output_CPI <- fit_JAGS(
  anomaly_detection_string,
  data_list_CPI,
  to_save,
  n_iter=10000
)

Once obtained both \(\{\delta_t\}_{t=1...N}\) sequences we can classify their values according to their expected values. Concretely, we classify a timestamp \(t\) such that \(\mathbb{E}[\delta_t] < 0.2\) as “high variance” (where strong oscillations may occur) and timestamps \(t\) with \(0.2 < \mathbb{E}[\delta_t] < 0.6\) as samples with “medium variance” (where we still might have an anomaly, but of weaker nature). All remaining samples can be considered stable and only subject to natural stochastic oscillations.

# Classifying anomalies for both series
anomalies_GDP <- classify_anomalies(anomaly_output_GDP, N_tot, title="GDP")

## Found 36 strong and 27 medium outliers out 305 total samples for GDP

anomalies_CPI <- classify_anomalies(anomaly_output_CPI, N_tot, title="CPI")

## Found 24 strong and 45 medium outliers out 305 total samples for CPI

# Displaying graphically the classification
delta_GDP <- anomaly_output_GDP$mean$d
delta_CPI <- anomaly_output_CPI$mean$d

display_classification(
  x_axis,
  delta_GDP,
  anomalies_GDP$high_var_indices,
  anomalies_GDP$med_var_indices,
  anomalies_GDP$low_var_indices,
  title="GDP"
)

display_classification(
  x_axis,
  delta_CPI,
  anomalies_CPI$high_var_indices,
  anomalies_CPI$med_var_indices,
  anomalies_CPI$low_var_indices,
  title="CPI"
)

Finally, we can plot both our time series displaying the points in time where high variances were observed: red vertical lines display timestamps where strong anomalies were detected, while orange dotted lines show medium anomalies.

# Zoom value: we only plot starting from this point on so as not to clog up the graph
start <- 150

visualise_anomalies(
  x_axis,
  series_GDP,
  anomalies_GDP$high_var_indices,
  anomalies_GDP$med_var_indices,
  title="GDP",
  start=start
)

visualise_anomalies(
  x_axis,
  series_CPI,
  anomalies_CPI$high_var_indices,
  anomalies_CPI$med_var_indices,
  title="CPI",
  start=start
)